Elaborato per Geospatial Data Analysis (MD2SL 2024) - prof.ssa Bocci

Il PIL degli Stati africani

Creazione dataset

rm(list=ls())

library(haven)

Warning: package ‘haven’ was built under R version 4.3.3

# lettura dati
setwd("C:/Users/dario/Documents/Projects/Master/Geospatial/")
stat <- read_dta("data/slave_trade_QJE.dta")

# rinonimo la colonna del PIL Log per capita GDP - from Maddison (2003)
colnames(stat)[colnames(stat) == "ln_maddison_pcgdp2000"] = "ln_pcgdp"
colnames(stat)[colnames(stat) == "ln_coastline_area"] = "ln_coast_area"
colnames(stat)[colnames(stat) == "ln_avg_gold_pop"] = "ln_gold_pop"
colnames(stat)[colnames(stat) == "ln_avg_oil_pop"] = "ln_oil_pop"
colnames(stat)[colnames(stat) == "ln_avg_all_diamonds_pop"] = "ln_diamonds_pop"
colnames(stat)[colnames(stat) == "atlantic_distance_minimum"] = "atlantic_dist_min"
colnames(stat)[colnames(stat) == "indian_distance_minimum"] = "indian_dist_min"
colnames(stat)[colnames(stat) == "saharan_distance_minimum"] = "saharan_dist_min"
colnames(stat)[colnames(stat) == "red_sea_distance_minimum"] = "red_sea_dist_min"


# rimuovo variabili non necessarie
delenda = c(#'abs_latitude', # Absolute latitude
            #'longitude', #Longitude
            'ln_export_pop', #Log total slave exports normalized by historic population
            'island_dum', #Indicator variable for small islands
            #'region_n', # Region indicator: North
            'region_s', # Region indicator: South
            'region_w', # Region indicator: West
            'region_e', # Region indicator: East
            'region_c', # Region indicator: Central
            'ln_pop_dens_1400', # Log population density in 1400
            'ethnic_fractionalization', # Ethnic fractionalization
            'state_dev', # State development
            'land_area') # Land area in millions of square kms
stat[,delenda] = NULL

library(sf)

Warning: package ‘sf’ was built under R version 4.3.3Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

# lettura shapefile
world = st_read('data/countries/50m/ne_50m_admin_0_countries.shp')

Reading layer `ne_50m_admin_0_countries' from data source 
  `C:\Users\dario\Documents\Projects\Master\Geospatial\data\countries\50m\ne_50m_admin_0_countries.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 242 features and 168 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
Geodetic CRS:  WGS 84

africa = world[which(world$CONTINENT=='Africa'),c('ISO_A3', 'NAME')]

par(mar=c(0,0,0,9))
coords <- st_coordinates(st_centroid(africa))

Warning: st_centroid assumes attributes are constant over geometries

plot(st_geometry(africa), border='gray', col='snow')
text(coords[, "X"], coords[, "Y"], labels=africa$ISO_A3, cex=0.5)

legend(x=st_bbox(africa)$xmax+1, y=st_bbox(africa)$ymax, 
       legend = paste(africa$ISO_A3, ' ', africa$NAME), cex=0.5, ncol=2, xpd=T)
par(mar=c(5, 4, 4, 2)+0.1)

Origine del dataset

Esiste un legame causale tra le deportazioni degli schiavi (ln_export_area) ed il PIL procapite (ln_pcgdp) degli odierni Stati africani?

par(mar=c(4.5,4,0.5,0.5))
plot(x=stat$ln_export_area, y=stat$ln_pcgdp, pch=16, bty='n', cex=0.7,
     xlab='ln_export_area', ylab='ln_pcgdp', col='firebrick')
text(x=stat$ln_export_area, y=stat$ln_pcgdp+0.1, labels=stat$isocode, cex=0.5)

grid()
par(mar=c(5,4,4,2)+0.1)

Rimozione isole

#rimuovo le isole
delenda = c('SYC', # Seychelles
            'STP', # Sao Tome & Principe
            'MUS', # Mauritius
            'MDG', #    Madagascar
            'CPV', #    Cape Verde Islands
            'COM'  #    Comoros
)
stat = stat[!(stat$isocode %in% delenda), ]
africa = africa[!(africa$ISO_A3 %in% delenda),]

Ricostruzione degi vecchi Stati africani

Rispetto a quando i dati sul PIL sono stati raccolti (fine anni ’90), la conformazione politica dell’Africa è cambiata; alcuni Stati si sono separati, per cui è necessario riunificare le loro mappe attuali:

• Sudan e Sudan del Sud riuniti nel Sudan
• Eritrea e Etiopia riuniti in Etiopia
• Somaliland e Somalia riunificati nella Somalia
• Sahara Occidentale e Marocco fusi nel Marocco

merge_countries <- function(sf_data, iso1, iso2, merged_iso, merged_name)
{
   country1 = sf_data[which(sf_data$ISO_A3 == iso1), ]
   country2 = sf_data[which(sf_data$ISO_A3 == iso2), ]
   merged_geometry = st_union(st_geometry(country1), st_geometry(country2))
   
   merged_data = data.frame(ISO_A3=merged_iso, NAME=merged_name, geometry=merged_geometry)
   merged_sf = st_sf(merged_data, crs=st_crs(sf_data))
   
   sf_data_filtered = sf_data[!(sf_data$ISO_A3 %in% c(iso1, iso2)), ]
   sf_data_merged = rbind(sf_data_filtered, merged_sf)
   return(sf_data_merged)
}


sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])

(africa[tofix_shape,c('NAME', 'ISO_A3')])

Simple feature collection with 5 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -17.09878 ymin: -13.45381 xmax: 48.93857 ymax: 27.65645
Geodetic CRS:  WGS 84
               NAME ISO_A3                       geometry
52         S. Sudan    SSD MULTIPOLYGON (((33.97607 4....
59       Somaliland    -99 MULTIPOLYGON (((48.93857 11...
105       W. Sahara    ESH MULTIPOLYGON (((-8.817773 2...
174         Eritrea    ERI MULTIPOLYGON (((36.52432 14...
192 Dem. Rep. Congo    COD MULTIPOLYGON (((30.75117 -8...

# Rinonimo Dem. Rep. Congo (COD) in Democratic Republic of Congo (ZAR)
africa$ISO_A3[which(africa$ISO_A3=='COD')] = 'ZAR'
africa$NAME[which(africa$ISO_A3=='COD')] = 'Democratic Republic of Congo'

# Eritrea -> Etiopia
africa = merge_countries(africa, 'ERI', 'ETH', 'ETH', 'Ethiope')
sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])

(africa[tofix_shape, c('NAME', 'ISO_A3')])

Simple feature collection with 3 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -17.09878 ymin: 3.490723 xmax: 48.93857 ymax: 27.65645
Geodetic CRS:  WGS 84
          NAME ISO_A3                       geometry
52    S. Sudan    SSD MULTIPOLYGON (((33.97607 4....
59  Somaliland    -99 MULTIPOLYGON (((48.93857 11...
105  W. Sahara    ESH MULTIPOLYGON (((-8.817773 2...

# W. Sahara -> Morocco
africa = merge_countries(africa, 'ESH', 'MAR', 'MAR', 'Morocco')
sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])

(africa[tofix_shape, c('NAME', 'ISO_A3')])

Simple feature collection with 2 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 24.14736 ymin: 3.490723 xmax: 48.93857 ymax: 12.2231
Geodetic CRS:  WGS 84
         NAME ISO_A3                       geometry
52   S. Sudan    SSD MULTIPOLYGON (((33.97607 4....
59 Somaliland    -99 MULTIPOLYGON (((48.93857 11...

# Somaliland -> Somalia
africa = merge_countries(africa, '-99', 'SOM', 'SOM', 'Somalia')
sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])

(africa[tofix_shape, c('NAME', 'ISO_A3')])

Simple feature collection with 1 feature and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 24.14736 ymin: 3.490723 xmax: 35.26836 ymax: 12.2231
Geodetic CRS:  WGS 84
       NAME ISO_A3                       geometry
52 S. Sudan    SSD MULTIPOLYGON (((33.97607 4....

# S. Sudan -> Sudan
africa = merge_countries(africa, 'SSD', 'SDN', 'SDN', 'Sudan')
sameCountries = intersect(africa$ISO_A3, stat$isocode)
tofix_stat = which(!stat$isocode %in% sameCountries)
tofix_shape = which(!africa$ISO_A3 %in% sameCountries)
(stat[tofix_stat, c('country','isocode')])

(africa[tofix_shape,c('NAME', 'ISO_A3')])

Simple feature collection with 0 features and 2 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
[1] NAME     ISO_A3   geometry
<0 rows> (or 0-length row.names)

Rimozione delle isolette di ogni Stato

Rimuoviamo le isole appartenenti ad alcuni Stati (in particolare un arcipelago molto a Sud sotto la nazionalità del Sudafrica)

### Rimozione isole
remove_islands <- function(sf, iso_a3)
{
   country = sf[which(sf$ISO_A3 == iso_a3), ]
   
   if (nrow(country) > 0)
   {
      
      geom_type <- st_geometry_type(country)[1] # Get the geometry type of the FIRST feature
      
      if (geom_type %in% c("MULTIPOLYGON"))
      {
         #country_polygons <- st_cast(country, "POLYGON")
         country_polygons <- suppressWarnings(st_cast(country, "POLYGON")) # Suppress warning here

         areas = st_area(country_polygons)
         largest_polygon = country_polygons[which.max(areas),]

         # Extract the geometry of the largest polygon
         country_geom = st_geometry(largest_polygon)[[1]]
         # Create a new sf object for the mainland
         country_mainland_sf = st_sf(data.frame(ISO_A3 = iso_a3, NAME = country$NAME[1]),
                                    geometry = st_sfc(country_geom, crs = st_crs(country)))

         sf_filtered = sf[!(sf$ISO_A3 %in% iso_a3), ]
         sf_final <- rbind(sf_filtered, country_mainland_sf)
         
         return(sf_final)
      }
   }
   
   return(sf)
}


# Rimozione delle isole di ogni stato
delenda = c('TUN', # Tunisia
            'TZA', # Tanzania
            'ZAF', # S. Africa
            'SLE', # Sierra Leone
            'MRT', # Mauritania
            'MWI', # Malawi
            'GNB', # Guinea-Bissau
            'GNQ', # Eq. Guinea
            'ETH') # Ethiope

par(mfrow = c(1, 2), mar=c(0,0,1,0))
for (c in delenda)
{
   plot(st_geometry(africa[which(africa$ISO_A3==c),'ISO_A3']), main=c, col='snow')
   africa = remove_islands(africa, c)
   plot(st_geometry(africa[which(africa$ISO_A3==c),'ISO_A3']), main=c, col='snow')
}

par(mfrow = c(1,1), mar=c(5, 4, 4, 2)+0.1)

Merge dei dataset

names(africa)[1:2] = names(stat[1:2])
stat$country = NULL #non mi serve più
ds = merge(africa, stat, by = "isocode")

coords <- st_coordinates(st_centroid(ds))

Warning: st_centroid assumes attributes are constant over geometries

par(mar=c(0,0,0,8))
plot(st_geometry(ds), border='gray', col='snow')
text(coords[, "X"], coords[, "Y"], labels=ds$isocode, cex=0.5)

legend(x=st_bbox(ds)$xmax+1, y=st_bbox(ds)$ymax, 
       legend = paste(ds$isocode, ' ', ds$country), cex=0.5, ncol=2, xpd=T)
par(mar=c(5, 4, 4, 2)+0.1)

Proiezione cartografica e conversione in SpatialPolygonsDataFrame

EPSG:2312 Garoua / UTM zone 33N (https://spatialreference.org/ref/epsg/2312/)

#utm <- st_crs("+proj=utm +zone=33 +north +datum=WGS84 +units=m +no_defs")
utm <- st_crs(2312)
sf_utm <- st_transform(ds, crs = utm)
st_crs(sf_utm)

Coordinate Reference System:
  User input: EPSG:2312 
  wkt:
PROJCRS["Garoua / UTM zone 33N",
    BASEGEOGCRS["Garoua",
        DATUM["Garoua",
            ELLIPSOID["Clarke 1880 (RGS)",6378249.145,293.465,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4197]],
    CONVERSION["UTM zone 33N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",15,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Cameroon - Garoua area."],
        BBOX[8.92,12.9,9.87,14.19]],
    ID["EPSG",2312]]

par(mfrow = c(1, 2), mar=c(0,0,2,0))
plot(st_geometry(ds), col = "lightcyan", main = "Original Projection")
plot(st_geometry(sf_utm), col = "lightyellow", main = "EPSG:2312 Projection")

par(mfrow = c(1, 1), mar=c(5, 4, 4, 2)+0.1)

# Lo salvo come shapefile per aprirlo in Geoda
# Save as Shapefile
st_write(sf_utm, "./data/slaveTrade.shp", driver="ESRI Shapefile", delete_layer=T) 

Warning: Field names abbreviated for ESRI Shapefile driver

Deleting layer `slaveTrade' using driver `ESRI Shapefile'
Writing layer `slaveTrade' to data source 
  `./data/slaveTrade.shp' using driver `ESRI Shapefile'
Writing 46 features with 29 fields and geometry type Unknown (any).

# Lo converto in SpatialPolygonDataFrame per ri-usare codice prof
slaveTrade = as_Spatial(sf_utm)

Descrizione del dataset

• isocode Country isocode
• country Country name
• ln_pcgdp Log per capita GDP - from Maddison (2003)
• ln_export_area Log total slave exports normalized by land area
• colony0 Colonizer indicator: not colonized
• colony1 Colonizer indicator: Britan
• colony2 Colonizer indicator: France
• colony3 Colonizer indicator: Portugal
• colony4 Colonizer indicator: Belgium
• colony5 Colonizer indicator: Spain
• colony6 Colonizer indicator: UN
• colony7 Colonizer indicator: Italy
• abs_latitude Absolute latitude
• longitude Longitude
• rain_min Min of monthly average rainfall (mm)
• humid_max Max of monthly afternoon avg humidity (%)
• low_temp Min of avg monthly low temp (C)
• ln_coastline_area Log (coastline/land_area)
• islam Percent Islamic
• legor_fr Legal origin indicator: French
• legor_uk Legal origin indicator: British
• region_n Region indicator: North Africa
• ln_avg_gold_pop Log gold production per capita
• ln_avg_oil_pop Log oil production per capita
• ln_avg_all_diamonds_pop Log diamond production per capita
• atlantic_distance_minimum Minimum Atlantic distance (000s of kms)
• indian_distance_minimum Minimum Indian distance (000s of kms)
• saharan_distance_minimum Minimum Saharan distance (000s of kms)
• red_sea_distance_minimum Minimum Red Sea distance (000s of kms)

library(sp)
library(spdep)

Warning: package ‘spdep’ was built under R version 4.3.3Loading required package: spData
Warning: package ‘spData’ was built under R version 4.3.3To access larger datasets in this package, install the spDataLarge package
with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`

Attaching package: ‘spData’

The following objects are masked _by_ ‘.GlobalEnv’:

    coords, world

library(spatialreg)

Warning: package ‘spatialreg’ was built under R version 4.3.3Loading required package: Matrix

Attaching package: ‘spatialreg’

The following objects are masked from ‘package:spdep’:

    get.ClusterOption, get.coresOption, get.mcOption, get.VerboseOption,
    get.ZeroPolicyOption, set.ClusterOption, set.coresOption, set.mcOption,
    set.VerboseOption, set.ZeroPolicyOption

library(RColorBrewer)
library(fields) # For the colorbar

Warning: package ‘fields’ was built under R version 4.3.3Loading required package: spam
Warning: package ‘spam’ was built under R version 4.3.3Spam version 2.11-0 (2024-10-03) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction 
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

Attaching package: ‘spam’

The following object is masked from ‘package:Matrix’:

    det

The following objects are masked from ‘package:base’:

    backsolve, forwardsolve

Loading required package: viridisLite
Warning: package ‘viridisLite’ was built under R version 4.3.3
Try help(fields) to get started.

library(gridExtra)

Warning: package ‘gridExtra’ was built under R version 4.3.3

matlab.like.hot <- function(n) 
{
   
   #my_gradient = c('gray25', 'red', 'yellow', 'floralwhite')
   my_gradient = c('gray25', 'red', 'yellow', 'cornsilk')
   #my_gradient = c('gray25', 'red', 'yellow', 'lightgoldenrodyellow')
   #my_gradient = c('gray25', 'red', 'yellow', 'lightyellow')
   #my_gradient = c('gray25', 'red', 'yellow', 'lemonchiffon')
   #my_gradient = c('gray25', 'red', 'yellow', 'snow')
   #my_gradient = c('gray25', 'red', 'yellow', 'ivory')
   return (colorRampPalette(my_gradient, space = "Lab")(n))
}

plotMapData <- function(shp, vname, min_v=NA, max_v=NA)
{
   n_colors <- 16
   # spectral_colors = brewer.pal(11, "Spectral")
   # spectral_r_palette <- colorRampPalette(rev(spectral_colors)) # Reverse for _r
   palette = matlab.like.hot(n_colors) #spectral_r_palette(n_colors)
   
   idv = which(names(shp)==vname)
   vals = shp[[idv]]
   
   if (is.na(min_v)) (min_v = min(vals))
   if (is.na(max_v)) (max_v = max(vals))
   if (min_v == max_v)
   {
      (min_v = floor(min(vals)))
      (max_v = ceiling(max(vals)))
   }
   
   norm_v = ((vals - min_v) / (max_v - min_v))
   norm_v[norm_v < 0] = 0
   norm_v[norm_v > 1] = 1
   pcol = palette[round(norm_v * (n_colors-1)) + 1]
   
   par(mar = c(0.0, 0.0, 2.0, 3))
   plot(shp, col = pcol,  lwd=0.5, main=vname, border='lightskyblue4')
   # Colorbar
   image.plot(legend.only = TRUE, zlim = c(min_v, max_v),
           col = palette,
           axis.args=list(cex.axis=0.8))
   par(mar=c(5,4,4,2) + 0.1)# Reset to default layout

   # spplot(shp, vname, main=vname, lwd=0.5,
   #        col = 'lightskyblue4',
   #     col.regions = palette,
   #     par.settings = list(axis.line = list(col = NA)))

   
}

showData <- function(shp, col_name, min_v=NA, max_v=NA)
{
   #par(mfrow = c(1, 2), mar = c(3, 4, 1.0, 0.5))
   par(mfrow = c(1, 2), mar = c(2, 4, 1.0, 0.5))
   
   h = hist(shp[[col_name]],main=col_name,xlab='', 
        col='lightskyblue3', border='lightskyblue4')#, breaks='FD',)
   grid(lty='solid', lwd=0.5, col='white', nx=NA, ny=NULL)
   
   plotMapData(shp, col_name)
   
   par(mfrow=c(1,1), mar=c(5,4,4,2)+0.1)
}

#--------------------------------------------------------------------



#ln_maddison_pcgdp2000
showData(slaveTrade,'ln_pcgdp')

#ln_export_area
showData(slaveTrade,'ln_export_area')

# rain_min
showData(slaveTrade,'rain_min')

# humid_max
showData(slaveTrade,'humid_max')

# low_temp
showData(slaveTrade,'low_temp')

# ln_coastline_area
showData(slaveTrade,'ln_coast_area')

# islam
showData(slaveTrade,'islam')

# ln_avg_gold_pop
showData(slaveTrade,'ln_gold_pop')

# ln_avg_oil_pop
showData(slaveTrade,'ln_oil_pop')

# ln_avg_all_diamonds_pop
showData(slaveTrade,'ln_diamonds_pop')

# atlantic_distance_minimum
showData(slaveTrade,'atlantic_dist_min')

# indian_distance_minimum
showData(slaveTrade,'indian_dist_min')

# saharan_distance_minimum
showData(slaveTrade,'saharan_dist_min')

# red_sea_distance_minimum
showData(slaveTrade,'red_sea_dist_min')

#..........................................................
colony_codes = c('--', 'GB', 'FR', 'PT', 'BE', 'ES', 'UN', 'IT')
#colony_names = names(slaveTrade)[5:12]
colony_names = c('colony0', 'colony1', 'colony2', 'colony3', 'colony4', 'colony5', 'colony6', 'colony7')
colony_id = max.col(stat[,colony_names])
slaveTrade$colony = as.factor(colony_codes[colony_id])

legor_codes = c('FR', 'GB', '--')
legor_names = c('legor_fr', 'legor_uk')
legor_id = max.col(stat[,legor_names])
slaveTrade$legor = legor_codes[legor_id]
slaveTrade$legor[rowSums(stat[,legor_names]) == 0] = '--'
slaveTrade$legor = as.factor(legor_codes[legor_id])

library(gridExtra)
#pL = 
spplot(slaveTrade, 'colony', main='colonizers',  lwd=0.5,
            col.regions = brewer.pal(n=length(colony_codes), name='Pastel2'), par.settings = list(axis.line = list(col = NA)))

# pR = 
spplot(slaveTrade, 'legor', main='legislative origin', lwd=0.5,
            col.regions = brewer.pal(n=length(legor_codes), name='Pastel1'), par.settings = list(axis.line = list(col = NA)))

#grid.arrange(pL, pR, ncol = 2)

Analisi di correlazione (variabili continue)

library(corrplot)

Warning: package ‘corrplot’ was built under R version 4.3.3corrplot 0.95 loaded

#(preds = names(slaveTrade)[c(3:4,13:18, 20:26)])
preds = c('ln_export_area', 'abs_latitude', 'longitude',
          'rain_min', 'humid_max', 'low_temp', 'ln_coast_area',
          'islam', 'ln_gold_pop', 'ln_oil_pop', 'ln_diamonds_pop',
          'atlantic_dist_min', 'indian_dist_min', 
          'saharan_dist_min', 'red_sea_dist_min')
corpredictor = cor(slaveTrade@data[,preds], method="pearson")

par(mar=c(0,0,0,0))
corrplot(corpredictor, type = "upper", col=rev(colorRampPalette(brewer.pal(n=11, name='RdBu'))(100)), tl.col = "black", tl.srt = 45, tl.cex=0.8)
par(mar=c(5,4,4,2)+0.1)

(corpredictor['abs_latitude', 'low_temp'])

[1] -0.7768648

(corpredictor['longitude', c('atlantic_dist_min', 'indian_dist_min', 'red_sea_dist_min')])

atlantic_dist_min   indian_dist_min  red_sea_dist_min 
        0.8145862        -0.7918568        -0.7828318 

(corpredictor['islam', 'saharan_dist_min'])

[1] -0.7350223

(corpredictor['atlantic_dist_min', 'red_sea_dist_min'])

[1] -0.8306802

Ispezione geospaziale dei dati

Determinazione dei vicini

slaveTrade.nb <- poly2nb(slaveTrade)
slaveTrade.lw <- nb2listw(slaveTrade.nb, style = "W")
slaveTrade.nb.lag = nblag(slaveTrade.nb, maxlag=5)


summary(slaveTrade.nb)

Neighbour list object:
Number of regions: 46 
Number of nonzero links: 202 
Percentage nonzero weights: 9.546314 
Average number of links: 4.391304 
Link number distribution:

1 2 3 4 5 6 7 8 9 
2 8 7 7 8 8 2 3 1 
2 least connected regions:
17 23 with 1 link
1 most connected region:
44 with 9 links

summary(slaveTrade.lw)

Characteristics of weights list object:
Neighbour list object:
Number of regions: 46 
Number of nonzero links: 202 
Percentage nonzero weights: 9.546314 
Average number of links: 4.391304 
Link number distribution:

1 2 3 4 5 6 7 8 9 
2 8 7 7 8 8 2 3 1 
2 least connected regions:
17 23 with 1 link
1 most connected region:
44 with 9 links

Weights style: W 
Weights constants summary:

cnb = coordinates(slaveTrade)
par(mfrow=c(1,1), mar=c(0,0,1,0))
plot(slaveTrade, border='gray', lwd=0.5, main='prox. lag 1')
plot(slaveTrade.nb.lag[[1]], cnb, add=T, col='red', pch=16, cex=0.5, lwd=0.5)

#
par(mfrow=c(2,2), mar=c(0,0,1,0))

#
plot(slaveTrade, border='gray', main='prox. lag 2', lwd=0.5)
plot(slaveTrade.nb.lag[[2]], cnb, add=T, col='red', pch=16, cex=0.3, lwd=0.5)

#
plot(slaveTrade, border='gray', main='prox. lag 3', lwd=0.5)
plot(slaveTrade.nb.lag[[3]], cnb, add=T, col='red', pch=16, cex=0.3, lwd=0.5)

#
plot(slaveTrade, border='gray', main='prox. lag 4', lwd=0.5)
plot(slaveTrade.nb.lag[[4]], cnb, add=T, col='red', pch=16, cex=0.3, lwd=0.5)

#
plot(slaveTrade, border='gray', main='prox. lag 5', lwd=0.5)
plot(slaveTrade.nb.lag[[5]], cnb, add=T, col='red', pch=16, cex=0.3, lwd=0.5)

par(mfrow=c(2,2), mar=c(5, 4, 4, 2)+0.1)

Tabella complessiva dei global Moran per tutte le variabili

library(pgirmess)

moran_list = list()
#cvars = names(slaveTrade)[c(3,4,13:26)]
cvars = c("ln_pcgdp", "ln_export_area", "rain_min", "humid_max", "low_temp", "ln_coast_area",
"islam", 
"ln_gold_pop","ln_oil_pop", "ln_diamonds_pop",
"atlantic_dist_min", "indian_dist_min", "saharan_dist_min", "red_sea_dist_min")

for (varname in cvars)
{
   vals = slaveTrade[[varname]]
   resM = moran.test(vals, slaveTrade.lw)
   resG = geary.test(vals, slaveTrade.lw)
   moran_list = rbind(moran_list, c(varname, resM$estimate[1], resM$p.value, resG$estimate[1], resG$p.value))
   
corD<-correlog(coordinates(slaveTrade), slaveTrade[[varname]], method="Moran")

par(mar=c(4, 3, 2, 0))
barplot(corD[,'coef'], names.arg=1:length(corD[,'coef']), yaxt = "n",
        main=paste('correlog. ', varname), xlab='lag', ylim=c(-1,1),
        col='lightskyblue3', border='lightskyblue4')
axis(2,  # 2 indicates the left axis (y-axis)
     at = seq(-1,1,0.2),  # Positions of the ticks
     labels = seq(-1,1,0.2))
#grid(lty='solid', lwd=0.5, col='white', nx=NA, ny=NULL)
abline(h = seq(-1,1,0.2), col = "white", lty = 'solid', lwd=0.5) # Add gridlines at y_ticks
par(mar=c(5, 4, 4, 2)+0.1)
print(varname)

}

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "ln_pcgdp"

[1] "ln_export_area"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "rain_min"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "humid_max"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "low_temp"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "ln_coast_area"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "islam"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "ln_gold_pop"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "ln_oil_pop"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "ln_diamonds_pop"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "atlantic_dist_min"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "indian_dist_min"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "saharan_dist_min"

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

[1] "red_sea_dist_min"

print(moran_list)

                          Moran I statistic                         
 [1,] "ln_pcgdp"          "0.423212218627437" "5.84693606250667e-06"
 [2,] "ln_export_area"    "0.305517058870406" "0.000707523106463913"
 [3,] "rain_min"          "0.164251212033616" "0.0270892478720973"  
 [4,] "humid_max"         "0.273673184682141" "0.00146370187499353" 
 [5,] "low_temp"          "0.569670215279815" "3.870356674188e-09"  
 [6,] "ln_coast_area"     "0.242017671557775" "0.00510666852032284" 
 [7,] "islam"             "0.705157163497059" "7.96725428967743e-13"
 [8,] "ln_gold_pop"       "0.174254358907525" "0.0284356869275728"  
 [9,] "ln_oil_pop"        "0.467275141106285" "7.95651993408283e-07"
[10,] "ln_diamonds_pop"   "0.313794535588015" "0.000427916545590119"
[11,] "atlantic_dist_min" "0.725556284996984" "6.49073307389046e-14"
[12,] "indian_dist_min"   "0.803117454548449" "3.44740790902477e-16"
[13,] "saharan_dist_min"  "0.87224731941177"  "6.07174170092653e-19"
[14,] "red_sea_dist_min"  "0.871485990898349" "1.06323009420501e-18"
      Geary C statistic                          
 [1,] "0.562737299594007"  "5.5455869476088e-05" 
 [2,] "0.551065826441476"  "1.47133692708474e-05"
 [3,] "0.799389686320364"  "0.0686756946949285"  
 [4,] "0.647335019269922"  "0.00217855925159111" 
 [5,] "0.326450740075291"  "2.62229760519281e-10"
 [6,] "0.683154733931286"  "0.00145453722187371" 
 [7,] "0.286658983545689"  "8.78501766139294e-12"
 [8,] "0.801654331091"     "0.0291452887728903"  
 [9,] "0.522042351663106"  "8.61703934180597e-06"
[10,] "0.674172952769783"  "0.00274798696239173" 
[11,] "0.242008049696664"  "3.78528596268506e-11"
[12,] "0.205462740496694"  "2.38830906084901e-13"
[13,] "0.0968316131722766" "1.04236233334834e-15"
[14,] "0.105249867609864"  "2.72806775219181e-16"

showMoran <- function(shf, lmii, tlab, roi)
{
   n_colors = 256
   my_gradient = c('deepskyblue3', 'snow', 'firebrick3')
   palette = colorRampPalette(my_gradient, space = "Lab")(n_colors)

   
   vals = lmii[,1]
   #min_v = min(vals); max_v = max(vals);
   ref_v = 1;#max(abs(vals))
   min_v = -ref_v; max_v = ref_v;
   norm_v = (vals-min_v)/(max_v-min_v)
   norm_v[norm_v < 0] = 0
   norm_v[norm_v > 1] = 1
   pcol = palette[round(norm_v * (n_colors-1)) + 1]
   
   #tcol = ifelse(lmii[,5]>0.05, "darkgrey", "black")
   tcol = ifelse(1:length(vals) %in% roi, 'black', 'darkgray')
   lcex = ifelse(1:length(vals) %in% roi, 0.7, 0.5)
   
   par(mar = c(0.0, 0.0, 1.0, 3))
   
   plot(shf, col=pcol, lwd=0.5, border='olivedrab', main=paste('local Moran: ', tlab))

   # Colorbar
   image.plot(legend.only = TRUE, zlim = c(min_v, max_v),
           col = palette, 
           axis.args=list(cex.axis=0.8))

   coords = coordinates(slaveTrade)
   text(coords[,1], coords[,2], labels=ds$isocode, cex=lcex, col=tcol)
   
   par(mar=c(5,4,4,2)+0.1)
   
}


showClusterMap <- function(ds, lmii, col_lab, significance_level=0.05)
{
   par(mar = c(0.0, 0.0, 1.0, 4))
   p_value = lmii[,5]
   Moran_I = lmii[,1]
   llist = c('HH',   'LL',     'HL',        'LH',        'NS')
   clist = c('firebrick', 'royalblue', 'indianred1', 'lightblue1', 'lightgray')
   lid <- ifelse(p_value < significance_level,
                        ifelse(Moran_I > 0,
                               ifelse(ds[[col_lab]] > mean(ds[[col_lab]]), 1, 2),
                               ifelse(ds[[col_lab]] > mean(ds[[col_lab]]), 3, 4)),
                        5)
   cluster = llist[lid]
   ccol = clist[lid]
   plot(ds, col=ccol, lwd=0.5, border='snow', main=paste('clusters: ', col_lab))
   legend("right", legend = llist, fill = clist, bty='n')
   par(mar=c(5,4,4,2)+0.1)
}

Variabile di outcome (ln_pcgdp)

Media

computeLocalAvg <- function(dsp, variable_name, wlist)
{
   local_averages <- numeric(nrow(dsp)) # Initialize a vector to store the results
   for (i in 1:nrow(dsp)) 
   {
      neighbors_i <- wlist$neighbours[[i]] #nblist[[i]] # Neighbors of area i
      weights_i <- wlist$weights[[i]] #wlist[[i]]
      
      if (length(neighbors_i) > 0) 
      { # Check if the area has neighbors (for islands or edge cases)
         neighbor_values <- dsp@data[neighbors_i, variable_name]  # values of the neighbors
         #print(neighbor_values)
         local_averages[i] <- sum(weights_i * neighbor_values)/sum(weights_i) # weighted average
      }
      else
      { # Handle cases with no neighbors (e.g., islands):
         local_averages[i] <- dsp@data[i, variable_name] 
         warning(paste("Area", i, "has no neighbors.")) # or print a message
      }
   }
   
   return(local_averages)
}


showMapByQuintiles <- function(dsp, variable_name, mydata)
{
   dsp_sf = st_as_sf(dsp)  # Convert to sf object
   dsp_sf$mydata <- mydata

   brks=quantile(mydata, seq(0, 1, 0.2), na.rm=T);
   brks[1] = floor(brks[1]*10)/10; 
   brks[length(brks)]=ceiling(brks[length(brks)]*10)/10
   spplot(as_Spatial(dsp_sf), 'mydata', main=paste('Quintile cuts: ', variable_name), lwd=0.5,
       at = brks,
       col = 'lightskyblue4',
       # col.regions = colorRampPalette(rev(brewer.pal(11, "Spectral")))(256),
       col.regions = matlab.like.hot(16),
       par.settings = list(axis.line = list(col = NA)))
}


# brks=quantile(slaveTrade$ln_pcgdp, seq(0, 1, 0.2), na.rm=T);
# brks[1] = floor(brks[1]*10)/10; 
# brks[length(brks)]=ceiling(brks[length(brks)]*10)/10
# spplot(slaveTrade, c('ln_pcgdp'), main='ln_pcgdp', lwd=0.5,
#        at = brks,
#        #col.regions = colorRampPalette(rev(brewer.pal(11, "Spectral")))(256),
#        col.regions = matlab.like.hot(16),
#        par.settings = list(axis.line = list(col = NA)))

ln_pcgdp_avg = computeLocalAvg(slaveTrade, 'ln_pcgdp', slaveTrade.lw)

dsp_dummy=slaveTrade; dsp_dummy$ln_pcgdp_avg = ln_pcgdp_avg
showData(dsp_dummy,'ln_pcgdp_avg')

showMapByQuintiles(slaveTrade, 'ln_pcgdp_avg', ln_pcgdp_avg)

moran.test(slaveTrade$ln_pcgdp, slaveTrade.lw)

    Moran I test under randomisation

data:  slaveTrade$ln_pcgdp  
weights: slaveTrade.lw    

Moran I statistic standard deviate = 4.3832, p-value = 5.847e-06
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.42321222       -0.02222222        0.01032717 

corD<-correlog(coordinates(slaveTrade), slaveTrade$ln_pcgdp, method="Moran")

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

corD = cbind(corD, corD[,1:2]); colnames(corD)[5:6] = c('my_Moran', 'my_pvalue')
for (k in 1:nrow(corD))
{
   if (k<=length(slaveTrade.nb.lag))
   {
      nb_lag = slaveTrade.nb.lag[[k]]
      lw_klag <- nb2listw(nb_lag, style = 'W')
      resk = moran.test(slaveTrade$ln_pcgdp, lw_klag)
      corD[k,5] = resk$estimate[1];
      corD[k,6] = resk$p.value;
   }
   else
   {
      corD[k,5] = NA
      corD[k,6] = NA
   }
}
(corD)

      dist.class         coef      p.value   n    my_Moran    my_pvalue
 [1,]     472379  0.581321114 1.430501e-05 102  0.42321222 5.846936e-06
 [2,]    1103686  0.214094436 8.420569e-03 214 -0.09532285 8.516463e-01
 [3,]    1734993 -0.145930298 8.950403e-01 218 -0.20962476 9.985113e-01
 [4,]    2366300 -0.244325064 9.950600e-01 264 -0.11027096 9.035166e-01
 [5,]    2997607 -0.122315364 8.783822e-01 262  0.13563682 2.773817e-02
 [6,]    3628914  0.155421053 1.812701e-02 256          NA           NA
 [7,]    4260221 -0.075661501 6.954287e-01 188          NA           NA
 [8,]    4891528 -0.285111943 9.964635e-01 208          NA           NA
 [9,]    5522835 -0.076877245 6.657695e-01 156          NA           NA
[10,]    6154142  0.130179157 9.734806e-02 126          NA           NA
[11,]    6785449  0.183756900 8.884203e-02  62          NA           NA
[12,]    7416756 -0.009985953 2.489103e-01  14          NA           NA

plot(x=1:nrow(corD), y=corD[,'coef'], col='black', xlim=c(0,nrow(corD)), bty='n'); grid()
points(x=1:nrow(corD), y=corD[,'my_Moran'], col='firebrick', pch=16, add=T)

Warning: "add" is not a graphical parameter

abline(a=0,b=0,lwd=1)

corD<-correlog(coordinates(slaveTrade), ln_pcgdp_avg, method="Moran")

Warning: neighbour object has 16 sub-graphsWarning: neighbour object has 4 sub-graphsWarning: neighbour object has 2 sub-graphsWarning: neighbour object has 6 sub-graphsWarning: neighbour object has 12 sub-graphsWarning: neighbour object has 19 sub-graphsWarning: neighbour object has 30 sub-graphsWarning: neighbour object has 41 sub-graphs

par(mar=c(4, 3, 2, 0))

barplot(corD[,'coef'], names.arg=1:length(corD[,'coef']), yaxt = "n",
        main=paste('correlog. ', 'ln_pcgdp_avg'), xlab='lag', ylim=c(-1,1),
        col='lightskyblue3', border='lightskyblue4')
axis(2,  # 2 indicates the left axis (y-axis)
     at = seq(-1,1,0.2),  # Positions of the ticks
     labels = seq(-1,1,0.2))

#grid(lty='solid', lwd=0.5, col='white', nx=NA, ny=NULL)
abline(h = seq(-1,1,0.2), col = "white", lty = 'solid', lwd=0.5) # Add gridlines at y_ticks
par(mar=c(5, 4, 4, 2)+0.1)

NA
NA

Local Moran

set.seed(42)
moran.plot(slaveTrade$ln_pcgdp, slaveTrade.lw,
           xlab='ln_pcgdp', asp=1, bty='n')

print(slaveTrade$isocode[c(5,19,23)])

[1] "BWA" "GNQ" "LSO"

lmii = localmoran(slaveTrade$ln_pcgdp, slaveTrade.lw)
# lmii = localmoran_perm(slaveTrade$ln_pcgdp, slaveTrade.lw, nsim=999, iseed=123456789)

showMoran(slaveTrade, lmii, 'ln_pcgdp', c(5,19,23))

showClusterMap(slaveTrade, lmii, 'ln_pcgdp', significance_level=0.05)

Trattamento (ln_export_area)

Media

ln_export_area_avg = computeLocalAvg(slaveTrade, 'ln_export_area', slaveTrade.lw)

dsp_dummy=slaveTrade; dsp_dummy$ln_export_area_avg = ln_export_area_avg
showData(dsp_dummy,'ln_export_area_avg')

showMapByQuintiles(slaveTrade, 'ln_export_area', ln_export_area_avg)

Local Moran

moran.test(slaveTrade$ln_export_area, slaveTrade.lw)

    Moran I test under randomisation

data:  slaveTrade$ln_export_area  
weights: slaveTrade.lw    

Moran I statistic standard deviate = 3.1916, p-value = 0.0007075
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.30551706       -0.02222222        0.01054508 

mr = moran.plot(slaveTrade$ln_export_area, slaveTrade.lw, 
           xlab='ln_export_area', asp=1, bty='n')

print(slaveTrade$isocode[c(40,43)])

[1] "TUN" "ZAF"

lmii = localmoran(slaveTrade$ln_export_area, slaveTrade.lw)
# lmii = localmoran_perm(slaveTrade$ln_export_area, slaveTrade.lw, nsim=999, iseed=123456789)
showMoran(slaveTrade, lmii, 'ln_export_area',  c(40,43))

showClusterMap(slaveTrade, lmii, 'ln_export_area')

Modelli di Regressione

fmla <- ln_pcgdp ~ ln_export_area + colony1 + colony2 +colony3 + colony4 + colony5 + colony6 + colony7 + rain_min + humid_max + low_temp + ln_coast_area + islam + legor_fr + region_n + ln_gold_pop + ln_oil_pop + ln_diamonds_pop

Modello lineare

modLM <- lm(fmla, data = slaveTrade)
summary(modLM)

Call:
lm(formula = fmla, data = slaveTrade)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.54419 -0.22749 -0.01176  0.17298  0.75527 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      6.409676   0.880992   7.276 7.96e-08 ***
ln_export_area  -0.090781   0.025850  -3.512  0.00158 ** 
colony1          0.947924   0.449902   2.107  0.04455 *  
colony2          1.043905   0.441921   2.362  0.02563 *  
colony3          0.867653   0.454028   1.911  0.06667 .  
colony4          0.257498   0.582354   0.442  0.66189    
colony5          1.559744   0.705872   2.210  0.03579 *  
colony6          1.570997   0.727731   2.159  0.03992 *  
colony7          0.469763   0.648013   0.725  0.47473    
rain_min         0.001658   0.007255   0.229  0.82096    
humid_max        0.016627   0.008935   1.861  0.07367 .  
low_temp        -0.053279   0.016356  -3.257  0.00303 ** 
ln_coast_area    0.099378   0.035313   2.814  0.00901 ** 
islam           -0.002329   0.002782  -0.837  0.40971    
legor_fr         0.010320   0.421089   0.025  0.98063    
region_n        -0.131477   0.379578  -0.346  0.73174    
ln_gold_pop      0.012940   0.013661   0.947  0.35193    
ln_oil_pop       0.077399   0.022307   3.470  0.00177 ** 
ln_diamonds_pop -0.021472   0.035204  -0.610  0.54701    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.404 on 27 degrees of freedom
Multiple R-squared:  0.8326,    Adjusted R-squared:  0.7209 
F-statistic: 7.459 on 18 and 27 DF,  p-value: 2.204e-06

Test di Moran sui residui

moran.test(residuals(modLM), slaveTrade.lw)

    Moran I test under randomisation

data:  residuals(modLM)  
weights: slaveTrade.lw    

Moran I statistic standard deviate = -0.36075, p-value = 0.6409
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      -0.05897598       -0.02222222        0.01037988 

lm.morantest(modLM, slaveTrade.lw)

    Global Moran I for regression residuals

data:  
model: lm(formula = fmla, data = slaveTrade)
weights: slaveTrade.lw

Moran I statistic standard deviate = 0.55695, p-value = 0.2888
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
    -0.058975984     -0.108068846      0.007769563 

Test di Rao

Anche se dal test di Moran risulta assenza di struttura spaziale sui residui, e pertanto quanto osservato su ln_pcgdp è già spiegato dalle covariate impiegate, applichiamo anche il test dei moltiplicatori di Lagrange (o test di Rao) per valutare la necessità di un modello a struttura spaziale, tramite lm.RStests().

La funzione lm.RStests() restituisce

• il test RSerr (LM-error, per SEM), che verifica l’autocorrelazione spaziale nel termine di errore;
• il test RSlag (LM-lag, per SLM), che verifica l’autocorrelazione spaziale nella variabile dipendente;
• le versioni robuste dei precedenti test;
• il test SARMA che verifica la dipendenza combinata tra lag spaziale ed errore.

#Rao's score (a.k.a Lagrange multiplier) diagnostics
lmtest = lm.RStests(modLM,listw = slaveTrade.lw, test="all" )
summary(lmtest)

    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence
data:  
model: lm(formula = fmla, data = slaveTrade)
test weights: slaveTrade.lw
 
         statistic parameter p.value
RSerr    0.3052046         1  0.5806
RSlag    0.0099526         1  0.9205
adjRSerr 0.4906222         1  0.4836
adjRSlag 0.1953702         1  0.6585
SARMA    0.5005748         2  0.7786

Gli elvati p-value forniti da tutti i test sembrano indicare che non esista nessuna dipendenza statisticamente significativa da una struttura spaziale.

Decidiamo comunque di applicare anche gli Spatial Error Model (SEM), Spatial Lag Model (SLM) e lo Spatial Durbin Model (SDM) per ulteriore conferma.

Spatial Error Model

library(spatialreg)
modSEM <- errorsarlm(fmla, data = slaveTrade, listw = slaveTrade.lw)
summary(modSEM)

Call:errorsarlm(formula = fmla, data = slaveTrade, listw = slaveTrade.lw)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.492646 -0.216874 -0.047111  0.218542  0.611678 

Type: error 
Coefficients: (asymptotic standard errors) 
                  Estimate Std. Error z value  Pr(>|z|)
(Intercept)      4.7163055  0.6755299  6.9816 2.918e-12
ln_export_area  -0.0599673  0.0176409 -3.3993 0.0006755
colony1          2.1542058  0.3550998  6.0665 1.307e-09
colony2          1.8236501  0.3169903  5.7530 8.767e-09
colony3          1.2430111  0.3353017  3.7071 0.0002096
colony4          1.1397081  0.3827486  2.9777 0.0029043
colony5          2.5782346  0.4802824  5.3682 7.954e-08
colony6          2.8080166  0.5338434  5.2600 1.441e-07
colony7          1.6344453  0.5735883  2.8495 0.0043787
rain_min         0.0082390  0.0047167  1.7468 0.0806799
humid_max        0.0230702  0.0067576  3.4140 0.0006402
low_temp        -0.0772995  0.0108697 -7.1115 1.148e-12
ln_coast_area    0.1478062  0.0235635  6.2727 3.549e-10
islam           -0.0022717  0.0017485 -1.2992 0.1938658
legor_fr         0.6507630  0.3589125  1.8132 0.0698084
region_n        -0.2911240  0.2332548 -1.2481 0.2119964
ln_gold_pop      0.0158570  0.0106232  1.4927 0.1355246
ln_oil_pop       0.0514046  0.0142822  3.5992 0.0003192
ln_diamonds_pop  0.0016549  0.0251344  0.0658 0.9475043

Lambda: -0.91746, LR test value: 2.344, p-value: 0.12577
Asymptotic standard error: 0.1741
    z-value: -5.2696, p-value: 1.3672e-07
Wald statistic: 27.769, p-value: 1.3672e-07

Log likelihood: -10.15064 for error model
ML residual variance (sigma squared): 0.075562, (sigma: 0.27489)
Number of observations: 46 
Number of parameters estimated: 21 
AIC: 62.301, (AIC for lm: 62.645)

modSEM$coefficients["lambda"]

<NA> 
  NA 

Test di Moran sui residui

moran.test(residuals(modSEM), slaveTrade.lw)

    Moran I test under randomisation

data:  residuals(modSEM)  
weights: slaveTrade.lw    

Moran I statistic standard deviate = -0.27332, p-value = 0.6077
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      -0.05021622       -0.02222222        0.01049021 

Spatial Lag Model

modSLM <- lagsarlm(fmla, data=slaveTrade, listw=slaveTrade.lw)
summary(modSLM)

Call:lagsarlm(formula = fmla, data = slaveTrade, listw = slaveTrade.lw)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.5410864 -0.2299147 -0.0046404  0.1681190  0.7515161 

Type: lag 
Coefficients: (asymptotic standard errors) 
                  Estimate Std. Error z value  Pr(>|z|)
(Intercept)      6.5317530  1.4034319  4.6541 3.254e-06
ln_export_area  -0.0912573  0.0213286 -4.2786 1.880e-05
colony1          0.9563989  0.3455217  2.7680 0.0056404
colony2          1.0558252  0.3397165  3.1080 0.0018838
colony3          0.8774812  0.3498822  2.5079 0.0121439
colony4          0.2609079  0.4499668  0.5798 0.5620239
colony5          1.5773265  0.5431796  2.9039 0.0036857
colony6          1.5721306  0.5627857  2.7935 0.0052144
colony7          0.4754505  0.4972550  0.9562 0.3389963
rain_min         0.0016561  0.0056141  0.2950 0.7680055
humid_max        0.0166923  0.0068453  2.4385 0.0147478
low_temp        -0.0542420  0.0135019 -4.0174 5.885e-05
ln_coast_area    0.1000657  0.0270624  3.6976 0.0002177
islam           -0.0023246  0.0021323 -1.0902 0.2756382
legor_fr         0.0114815  0.3240160  0.0354 0.9717329
region_n        -0.1394538  0.2951216 -0.4725 0.6365484
ln_gold_pop      0.0129179  0.0104645  1.2345 0.2170335
ln_oil_pop       0.0780160  0.0189756  4.1114 3.933e-05
ln_diamonds_pop -0.0208386  0.0272283 -0.7653 0.4440767

Rho: -0.016867, LR test value: 0.011338, p-value: 0.9152
Asymptotic standard error: 0.14846
    z-value: -0.11361, p-value: 0.90955
Wald statistic: 0.012908, p-value: 0.90955

Log likelihood: -11.31695 for lag model
ML residual variance (sigma squared): 0.095761, (sigma: 0.30945)
Number of observations: 46 
Number of parameters estimated: 21 
AIC: 64.634, (AIC for lm: 62.645)
LM test for residual autocorrelation
test value: 0.55075, p-value: 0.45801

Test di Moran sui residui

moran.test(residuals(modSLM), slaveTrade.lw)

    Moran I test under randomisation

data:  residuals(modSLM)  
weights: slaveTrade.lw    

Moran I statistic standard deviate = -0.31678, p-value = 0.6243
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      -0.05450516       -0.02222222        0.01038525 

Spatial Durbin Model

modSDM <- lagsarlm(fmla, data=slaveTrade, listw=slaveTrade.lw, type='mixed')
summary(modSDM)

Call:lagsarlm(formula = fmla, data = slaveTrade, listw = slaveTrade.lw, 
    type = "mixed")

Residuals:
       Min         1Q     Median         3Q        Max 
-0.3454587 -0.1044586 -0.0089223  0.0780763  0.2933857 

Type: mixed 
Coefficients: (asymptotic standard errors) 
                       Estimate  Std. Error z value  Pr(>|z|)
(Intercept)          4.02561036  3.09324858  1.3014 0.1931153
ln_export_area      -0.15226732  0.02643228 -5.7607 8.379e-09
colony1              1.64439400  0.42314240  3.8861 0.0001018
colony2              1.74069948  0.35216696  4.9428 7.700e-07
colony3              1.03617073  0.36543845  2.8354 0.0045766
colony4              0.20900672  0.43540879  0.4800 0.6312102
colony5              1.65660160  0.38651872  4.2860 1.820e-05
colony6              2.51255171  0.53000474  4.7406 2.131e-06
colony7              1.60481168  0.60658073  2.6457 0.0081530
rain_min            -0.02337578  0.00568558 -4.1114 3.932e-05
humid_max            0.03772204  0.00820499  4.5975 4.277e-06
low_temp            -0.02299780  0.01526461 -1.5066 0.1319108
ln_coast_area        0.15205182  0.02752320  5.5245 3.304e-08
islam               -0.00514750  0.00189003 -2.7235 0.0064594
legor_fr            -0.06708916  0.45885552 -0.1462 0.8837558
region_n            -0.76731534  0.24208578 -3.1696 0.0015265
ln_gold_pop          0.04240323  0.01125507  3.7675 0.0001649
ln_oil_pop           0.13613948  0.02163123  6.2937 3.101e-10
ln_diamonds_pop     -0.06507716  0.02402345 -2.7089 0.0067506
lag.ln_export_area   0.05496841  0.04638199  1.1851 0.2359684
lag.colony1          5.14461052  1.10174121  4.6695 3.019e-06
lag.colony2          4.22291041  0.93596123  4.5118 6.427e-06
lag.colony3          2.10441422  0.91541635  2.2989 0.0215129
lag.colony4          3.44869799  0.85956982  4.0121 6.018e-05
lag.colony5          6.41462091  1.37376814  4.6694 3.021e-06
lag.colony6          1.53701918  1.58396026  0.9704 0.3318647
lag.colony7          7.93422220  1.50120046  5.2853 1.255e-07
lag.rain_min         0.00019509  0.01326682  0.0147 0.9882675
lag.humid_max       -0.00484986  0.02456831 -0.1974 0.8435121
lag.low_temp        -0.06941803  0.03564798 -1.9473 0.0514964
lag.ln_coast_area    0.12894507  0.06492060  1.9862 0.0470115
lag.islam           -0.00049062  0.00427151 -0.1149 0.9085565
lag.legor_fr         1.04046912  1.22216067  0.8513 0.3945828
lag.region_n         0.00018129  0.70397041  0.0003 0.9997945
lag.ln_gold_pop      0.10107588  0.02706225  3.7349 0.0001878
lag.ln_oil_pop      -0.16916945  0.04853784 -3.4853 0.0004916
lag.ln_diamonds_pop  0.01798481  0.06894153  0.2609 0.7941924

Rho: -0.48882, LR test value: 3.8581, p-value: 0.049506
Asymptotic standard error: 0.20297
    z-value: -2.4084, p-value: 0.016023
Wald statistic: 5.8004, p-value: 0.016023

Log likelihood: 20.26546 for mixed model
ML residual variance (sigma squared): 0.02305, (sigma: 0.15182)
Number of observations: 46 
Number of parameters estimated: 39 
AIC: 37.469, (AIC for lm: 39.327)
LM test for residual autocorrelation
test value: 11.48, p-value: 0.0007035

Test di Moran sui residui

moran.test(residuals(modSDM), slaveTrade.lw)

    Moran I test under randomisation

data:  residuals(modSDM)  
weights: slaveTrade.lw    

Moran I statistic standard deviate = -0.81573, p-value = 0.7927
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      -0.10527475       -0.02222222        0.01036613 

Confronto degli AIC

cmpMat = c()
cmpMat[1] = AIC(modLM)
cmpMat[2] = AIC(modSEM)
cmpMat[3] = AIC(modSLM)
cmpMat[4] = AIC(modSDM)

names(cmpMat) = c('LM', 'SEM', 'SLM', 'SDM')
(cmpMat)

      LM      SEM      SLM      SDM 
62.64525 62.30128 64.63391 37.46908 

Conclusioni

• Tutti i modelli testati riportano un indice di Moran sui residui statisticamente non diverso da 0
• Ne consegue che in tutti i casi le covariate riescono a spiegare la struttura spaziale rilevata dall’indice di Moran sulla variabile dipendente %l’indice di Moran osservato ed inizialmente imputato ad lieve trend
• Il modello con AIC più basso è lo Spatial Durbin Model (SDM) Tuttavia, applicando il rasoio di Occam, possiamo optare per il modello più semplice di tutti: quello lineare (LM)